import math

# based on OGP Guidance notes 7-2 formula
def stereo_scale(latF, e):
    if latF < 0: # assume the south-pole case
        tF = math.tan(math.pi/4. + latF/2.) / (((1 + e*math.sin(lat))/(1-e*math.sin(lat)))**(e/2.))
    else:        # assume the north-pole case
        tF = math.tan(math.pi/4. - latF/2.) * (((1 + e*math.sin(lat))/(1-e*math.sin(lat)))**(e/2.))
    #print 'tF =', tF
    mF = math.cos(lat)/((1 - (e**2)*(math.sin(lat)**2))**0.5)
    #print 'mF =', mF
    return mF * ((((1 + e)**(1 + e)) * ((1-e)**(1-e)))**0.5) / (2.*tF)

def deg2rad(deg):
    return deg * math.pi / 180.

if __name__ == '__main__':
    lat = deg2rad(71)
    e = 0.081819191  # e for WGS84
    print stereo_scale(lat, e)
